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ABSTRACT 

Aims. We present a comparison between independent computer codes, modeling the physics and chemistry of interstellar photon 
dominated regions (PDRs). Our goal was to understand the mutual differences in the PDR codes and their effects on the physical and 
chemical structure of the model clouds, and to converge the output of different codes to a common solution. 

Methods. A number of benchmark models have been created, covering low and high gas densities n = 10 3 , 10 5 5 cirr 3 and far 
ultraviolet intensities^ = 10, 10 5 in units of the Draine field (FUV: 6 < hv < 13.6 eV). The benchmark models were computed in 
two ways: one set assuming constant temperatures, thus testing the consistency of the chemical network and photo-processes, and a 
second set determining the temperature self consistently by solving the thermal balance, thus testing the modeling of the heating and 
cooling mechanisms accounting for the detailed energy balance throughout the clouds. 

Results. We investigated the impact of PDR geometry and agreed on the comparison of results from spherical and plane-parallel 
PDR models. We identified a number of key processes governing the chemical network which have been treated differently in the 
various codes such as the effect of PAHs on the electron density or the temperature dependence of the dissociation of CO by cosmic 
ray induced secondary photons, and defined a proper common treatment. We established a comprehensive set of reference models 
for ongoing and future PDR model bench-marking and were able to increase the agreement in model predictions for all benchmark 
models significantly. Nevertheless, the remaining spread in the computed observables such as the atomic fine-structure line intensities 
serves as a warning that there is still a considerable uncertainty when interpreting astronomical data with our models. 

Key words. ISM: abundances - Astrochemistry - ISM: clouds - ISM: general - Radiative Transfer - Methods: numerical 



1. Introduction 

Interstellar photon dominated regions or photodissociation re- 
gions (PDRs) play an important role in modern astrophysics as 
they are responsible for many emission characteristics of the 
ISM, and dominate the infrared and sub-millimetre spectra of 
star formation regions and galaxies as a whole. Theoretical mod- 
els addressing the structure of PDRs have been available for 
approximately 30 years and have evolved into advanced com- 
puter codes accounting for a growing number of physical ef- 
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fects with increasing accuracy. These codes have been devel- 
oped with different goals in mind: some are geared to efficiently 
model a particular type of region, e.g. HII regions, protoplan- 
etary disks, planetary nebulae, diffuse clouds, etc.; others em- 
phasize a strict handling of the micro-physical processes in full 
detail (e.g. wavelength dependent absorption), but at the cost of 
increased computing time. Yet others aim at efficient and rapid 
calculation of large model grids for comparison with observa- 
tional data, which comes at the cost of pragmatic approxima- 
tions using effective rates rather than detailed treatment. As a 
result, the different models have focused on the detailed simula- 
tion of particular processes determining the structure in the main 
regions of interest while using only rough approximations for 
other processes. The model setups vary strongly among different 
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model codes. This includes the assumed model geometry, their 
physical and chemical structure, the choice of free parameters, 
and other details. Consequently it is not always straightforward 
to directly compare the results from different PDR codes. Taking 
into account that there are multiple ways of implementing physi- 
cal effects in numerical codes, it is obvious that the model output 
of different PDR codes can differ from each other. As a result, 
significant variations in the physical and chemical PDR struc- 
ture predicted by the various PDR codes can occur. This diver- 
gence would prevent a unique interpretation of observed data 
in terms of the parameters of the observed clouds. Several new 
facilities such as Herschel, SOFIA, APEX, ALMA, and others 
will become available over the next years and will deliver many 
high quality observations of line and dust continuum emission 
in the sub-millimeter and FIR wavelength regime. Many impor- 
tant PDR tracers emit in this range ([CII] (158/xm), [OI] (63 and 
146 fxm), [CI] (370 and 610 /zm), CO (650, 520, 57.8 fim), 
H2O, etc.). In order to reliably analyze these data we need a set 
of high quality tools, including PDR models that are well un- 
derstood and properly debugged. As an important preparatory 
step toward these missions an international cooperation between 
many PDR model groups was initiated. The goals of this PDR- 
benchmarking were: 

- to understand the differences in the different code results 

- to obtain (as much as possible) the same model output with 
every PDR code when using the same input 

- to agree on the correct handling of important processes 

- to identify the specific limits of applicability of the available 
codes 

To this end, a PDR-benchmarking workshop was held at the 
Lorentz Center in Leiden, Netherlands in 2004 to jointly work 
on these topics Q. In this paper we present the results from this 
workshop and the results originating from the follow-up activi- 
ties. A related w orkshop to test line radiativ e transfer codes was 
held in 1999 (see lvan Zadelhoff et al.Ll2002h . 

It is not the purpose of the benchmarking to present a pre- 
ferred solution or a preferred code. PDRs are found in a large 
variety of objects and under very different conditions. To this 
end, it was neither possible nor desirable to develop a generic 
PDR code, able to model every possible PDR. Furthermore, the 
benchmarking is not meant to model any 'real' astronomical ob- 
ject. The main purpose of this study is technical not physical. 
This is also reflected in the choice of the adopted incomplete 
chemical reaction network (see §|4j. 

In §|2] we briefly introduce the physics involved in PDRs, in 
§|3] we introduce some key features in PDR modeling. §|4] de- 
scribes the setup of the benchmark calculations and § |5]presents 
the results for a selection of benchmark calculations and gives a 
short review over the participating codes. In § [6] we discuss the 
results and summarize the lessons learned from the benchmark 
effort. A tabular overview of the individual code characteristics 
is given in the Appendix. 



2. The Physics of PDRs 

PDRs are traditionally defined as regions where H2 -non-ionizing 
far-ultraviolet photons from stellar sources control the gas heat- 
ing and chemistry. Any ionizing radiation is assumed to be ab- 
sorbed in the narrow ionization fronts located between adjacent 



HII regions and the PDRsE In PDRs the gas is heated by the 
far-ultraviolet radiation (FUV, 6 < hv < 13.6 eV, from the 
ambient UV field and from hot stars) and cooled via the emis- 
sion of spectral line radiation of atomic and molecular species 
and continuum emission by dust (Hollenbach & Tielens 1999, 
Sternberg 2004). The FUV photons heat the gas by means of 
photoelectric emission from grain surfaces and polycyclic aro- 
matic hydrocarbons (PAHs) and by collisional de-excitation of 
vibrationally excited H2 molecules. Additional contributions to 
the total gas heating comes from IT formation, dissociation 
of H2, dust-gas collisions in case of dust temperatures exceed- 
ing the gas temperature, cosmic ray heating, turbulence, and 
from chemical heating. At low visual extinction Ay into the 
cloud/PDR the gas is cooled by emission of atomic fine-structure 
lines, mainly [OI] 63/mi and [CII] 158/mi. At larger depths, mil- 
limeter, sub-millimeter and far-infrared molecular rotational-line 
cooling (CO, OH, H2, H2O) becomes important, and a correct 
treatment of the radiative transfer in the line cooling is critical. 
The balance between heating and cooling determines the local 
gas temperature. The local FUV intensity also influences the 
chemical structure, i.e. the abundance of the individual chemi- 
cal constituents of the gas. The surface of PDRs is mainly dom- 
inated by reactions induced by UV photons, especially the ion- 
ization and dissociation of atoms and molecules. With diminish- 
ing FUV intensity at higher optical depths more complex species 
may be formed without being radiatively destroyed immediately. 
Thus the overall structure of a PDR is the result of a complex in- 
terplay between radiative transfer, energy balance, and chemical 
reactions. 

3. Modeling of PDRs 

The history of PDR modeling dates back to the early 1 970' s 
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Black & Dalgarno, 1977) with steady state models for the 
transitions from H to H2 and from C + to CO. In the fol- 
lowing years a number of models, addressing the chemi- 
cal and thermal structure of clouds subject t o an incident 



flux o f FUV photons have been d e veloped dde Jong et al 



1980; iTielens & Hollenbachi 119851: [van Dishoeck & Black . 



1 19881: ISternberg & Dalgarn oi Il989t iHoflenbach et al.L 11991 
Le Bourlot et all 119931: IStorzer et all 1 19961) Additionally, a 
number of models, focusing on certain aspects of PDR 
physics and chemistry were developed, e.g. models ac- 
counting for time-dependent chemical ne tworks, models of 
clumped media, and turbulent PDR models (Hill & H ollenbach, 
1 19781: IWagenblast & Hartquistl [19881: Ide Boisanger et all Il992[ 
bertoldi_&prainel 119961: Lee et al.L 119961: [Hegmann & Kegel 
[l996 | ISpaaiisl fl996t iNeiad & Wagenbla"sl fl999t iRollig et all 
BOOa iBell et all 120051) . Standard PDR models generally do not 
account for dynamical properties of gas but there are some stud- 
ies that consider the advection problem ra t her th an the steady 
state approach (e. g. IStorzer & Hollenbach, [ 19981) . For a more 
detailed review see Hollenba ch & Tielensl (Q999). 

In order to numerically model a PDR it is necessary to com- 
pute all local properties of a cloud such as the relative abun- 
dances of the gas constituents together with their level popula- 
tions, temperature of gas and dust, gas pressure, composition of 
dust/PAHs, and many more. This local treatment is complicated 
by the radiation field which couples remote parts of the cloud. 
The local mean radiation field, which is responsible for photo- 
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chemical reactions, gas/dust heating, and excitation of molecules 
depends on the position in the cloud and the (wavelength depen- 
dent) absorption along the lines of sight toward this position. 
This non-local coupling makes numerical PDR calculations a 
CPU time consuming task. 

PDR modelers and observers approach the PDRs from oppo- 
site sides: PDR models start by calculating the local properties 
of the clouds such as the local CO density and the correspond- 
ing gas temperature and use these local properties to infer the 
expected global properties of the cloud like total emergent emis- 
sivities or fluxes and column densities. The observer on the other 
hand starts by observing global features of a source and tries to 
infer the local properties from that. The connection between lo- 
cal and global properties is complex and not necessarily unam- 
biguous. Large variations e.g. in the CO density at the surface 
of the cloud may hardly affect the overall CO column density 
due to the dominance of the central part of the cloud with a high 
density. If one is interested in the total column density it does 
not matter whether different codes produce a different surface 
CO density. For the interpretation of high-J CO emission lines, 
however, different CO densities in the outer cloud layers make a 
huge difference since high temperatures are required to produce 
high-J CO fluxes. Thus, if different PDR model codes deviate in 
their predicted cloud structures, this may affect the interpretation 
of observations and may prevent inference of the 'true' structure 
behind the observed data. To this end it is very important to un- 
derstand the origin of present differences in PDR model calcu- 
lations. Otherwise it is impossible to rule out alternative inter- 
pretations. The ideal situation, from the modelers point of view, 
would be a complete knowledge of the true local structure of 
a real cloud and their global observable properties. This would 
easily allow us to calibrate PDR models. Since this case is unob- 
tainable, we take one step back and apply a different approach: 
If all PDR model codes use exactly the same input and the same 
model assumptions they should produce the same predictions. 

Because of the close interaction between chemical and ther- 
mal balance and radiative transfer, PDR codes typically iter- 
ate through the following computation steps: 1) solve the local 
chemical balance to determine local densities, 2) solve the local 
energy balance to estimate the local physical properties like tem- 
peratures, pressures, and level populations, 3) solve the radiative 
transfer, 4) for finite models it is necessary to successively it- 
erate steps l)-3). Each step requires a variety of assumptions 
and simplifications. Each of these aspects c an be investigated to 
great d etail and complexity (see for example lvan Zadelhoff et all 
(2002) for a discussion of NLTE radiative transfer methods), but 
the explicit aim of the PDR comparison workshop was to under- 
stand the interaction of all computation steps mentioned above. 
Even so it was necessary to considerably reduce the model com- 
plexity in order to disentangle cause and effect. 

3. 1 . Description of Sensitivities and Pitfalls 

Several aspects of PDR modeling have shown the need for de- 
tailed discussion, easily resulting in misleading conclusions if 
not treated properly: 

3.1 .1 . Model Geometry 

Two common geometrical setups of model PDRs are shown in 
Figure Q] Most PDR models feature a plane-parallel geometry, 
illuminated either from one side or from both sides. This ge- 
ometry naturally suggests a directed illumination, perpendicu- 
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Fig. 1. Common geometrical setups of a model PDR. The sur- 
face of any plane-parallel or spherical cloud is illuminated either 
a) uni-directional or b) isotropically. 



directed vs. isotropic attenuation 
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Fig. 2. Comparison of attenuation of the mean intensity for the 
case of an uni-directional and isotropically illuminated medium. 
The solid line gives the attenuation due to uni-directional il- 
lumination, while the dashed line gives the attenuation for an 
isotropic FUV radiation where r means the optical depth per- 
pendicular to the surface of the cloud. 



lar to the cloud surface. This simplifies the radiative transfer 
problem significantly, since it is sufficient to account for just 
one line of sight, if w e ignore scattering out of the line of sight 
dFlannerv et al.L Il980t) . Since most plane-parallel PDR models 
are infinite perpendicular to the cloud depth z it is also straight- 
forward to account for an isotropic FUV irradiation within the 
pure 1-D formalism. For a spherical geometry one can exploit 
the model symmetry only for a FUV field isotropically imping- 
ing onto the cloud. In finite plane-parallel and spherical models 
iterations over the depth/radial structure are mandatory because 
radiation is coming from multiple directions, passing through 
cloud elements for which the physical and chemical structure 
and hence opacities have not been calculated in the same itera- 
tion step. To account for this 'backside' illumination it is essen- 
tial to iterate on the radiation field. 
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The most important quantity describing the radiation field 
in PDR models is the local mean intensity (or alternatively the 
energy density) as given by: 



= - f 



, c/Q [erg cm 2 s 1 



Hz 



sr 1 ] 



(l) 



with the specific intensity I v being averaged over the solid an- 
gle £1 No te that when re ferring to the amb ient FUV in units of 
Draine X (iDrainel 1 19781) or Habing Go (iHabingl 1 19681) fields, 
these are always given as averaged over An. If we place a model 
cloud of sufficient optical thickness in such an average FUV 
field, the resulting local mean intensity at the cloud edge is half 
the value of that without the cloud. 

The choice between directed and isotropic FUV fields 
directly influences the attenuation due to dust. In the uni- 
directional case the FUV intensity along the line of sight is at- 
tenuated according to exp(-r v ), where t v is the optical depth of 
the dust at frequency v. For pure absorption the radiative transfer 
equation becomes: 



dl v (jx, x) 
dx 



= ~Ky /y(jU, X) 



(2) 



with the cosine of the radiation direction /j. = cos 0, the cloud 
depth x, and the absorption coefficient k v , with the simple so- 
lution Jy/J v fl = exp(-r vj u) for a semi-infinite cloud. For the 
isotropic case, 7 v ,o(y") = Jv,o = const., integration of Eq.[2]leads 
to the second order exponential integral: 



J V /Jyfi = E 2 (Ty) 



Jo 



exp(-T yJ u) 



dfi 



(3) 



As seen in Figure [2] the attenuation with depth in the isotropic 
case is significantly different from the uni-directional case. A 
common way to describe the depth dependence of a particular 
quantity in PDRs is to plot it against Ay, which is a direct mea- 
sure of the traversed column of attenuating material. In order 
to compare the uni-directional and the isotropic case it is nec- 
essary to rescale them to the same axis. It is possible to define 
an effective Ay, e ff = -\n[E2(Ay k)]/k with k = ruv/Ay in the 
isotropic case, where Ay is the attenuation perpendicular to the 
surface and UV is in the range 6 < hv < 13.6. In this paper 
all results from spherical models are scaled to Ay, e ff . Figure [3] 
demonstrates the importance of scaling results to an appropri- 
ate Ay scale. It shows the local H2 photo-dissociation rate for 
two different FUV illumination geometries. The solid line rep- 
resents a standard uni-directional illumination perpendicular to 
the cloud surface as given in many standard plane-parallel PDR 
codes. The dashed line is the result from an isotropic illumina- 
tion plotted against the standard 'perpendicular' Ay. The offset 
to the uni-directional case is significant. After rescaling to an ap- 
propriate Ay eif both model results are in good agreement. Please 
note, that in general it is not possible to achieve perfect agree- 
ment as there is a spectrum involved with a spread of k values 
across the UV. 

The attenuation of FUV radiation is additionally compli- 
cated if we account for dus t scattering. Fo r a fu ll treatment 
by Legendre polynomials see iFlannerv et all dl980h . In case of 
small scattering angles g = (cos 0) « 1 the scattering can be ap- 
proximated by an effective forward attenuation t(1 - to), where 
oj is the scattering albedo. Thus, more material is needed to ob- 
tain the same attenuation as in the case without scattering. Hence 
a proper scaling of Ay is necessary. In case of clumped gas this 
becomes even more complex. The presence of stochastic den- 
sity fluctuations leads to a substantial reduction of the effective 
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Fig. 3. H2 photo-dissociation rates resulting from uni-directional 
FUV illuminated clouds compared to an isotropic illumination. 
The results from isotropic models are plotted vs. the perpendic- 
ular Ay and vs. Ay e g. 



optical depth as demonstrated by Hegmann & Kegell d2003l) . All 
this has to be considered when calculating the photodissocia- 
tion and photoionization rates, when the attenutation with depth 
is represented by simple exponential form s, exp(-£,Ay) (e.g. 
van Dishoeck, ll988t iRoberge et all Il99ll) . where the factor h 
accounts for the wavelength dependence of the photoprocess Q 

3.1.2. Chemistry 

PDR chemistry has be e n addr e ssed in detail by many a uthors 



dTielens & Hollenbach, 1985t Ivan Dishoeck & Blackl 



Hollenbach etaTl [1991 



1993; 



Jansen et al.i 
— it — 



. IFuente et al.L 1 19931 
1995L 



Le Bourlot et al., 
Sternberg & Dalgarnol 1 1 9951: 



Lee et all Il996t Ba kes & Tielens, 1998; Walmsl ey et al.Ll l99S 
Sav age & ZiurvsL l2004t iTevssier et all l2004t IFuente et all 
2005; iMeiierink & Spaansl [2005). These authors discuss nu- 
merous aspects of PDR chemistry in great detail and give a 
comprehensive overview of the field. Here we repeat some 
crucial points in the chemistry of PDRs in order to motivate the 
benchmark standardization and to prepare the discussion of the 
benchmark result. 

In PDRs photoprocesses are very important due to the 
high FUV intensity, as well as reactions with abundant hy- 
drogen atoms. The formation and destruction of H2, heav- 
ily influenced by the FUV field, is of major importance for 
the chemistry in PDRs. H2 forms on grain surfaces, a pro- 
cess which cr ucially depends on the temperatures of the gas 
and th e grains (Hollenbach & Sal peteilll97ll ; ICazaux & Tie lens. 
2004), which themselves depend on the local cooling and heat- 
ing, governed by the FUV. The photo-dissociation of H2 is a 
line absorption process and, thus is subject to effective shield- 
ing ([van Disho eck & Blackl [1988). This leads to a sharp transi- 
tion from atomic to molecular hydrogen once the H2 absorption 
lines are optically thick. The photo-dissociation of CO is also 
a line absorption process, additionally complicated by the fact 
that the broad H2 absorption lines overlap with CO absorption 
lines. Similar to H2 this leads to a transition from atomic carbon 



3 In this context the term photoprocess refers to either photodissoci- 
ation or photoionization. 
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to CO. For Ay < 1 carbon is predominantly present in ionized 
form. For an assumed FUV field of^f = 1, CO is formed at about 
Ay » 2. This results in the typical PDR stratification of H/ H2 
and C + / CI CO. The depth of this transition zone depends on the 
physical parameters but also on the contents of the chemical net- 
work: for example the inclusion of PAHs into the chemical bal- 
ance calcuhiti£ns_slr^s to sm aller Av, e ff 
(e.g. lLepp & Dalgarnol [19881: iBakes & Tielenslll998l) . 

The solution of the chemical network itself covers the de- 
struction and formation reactions of all chemical species consid- 
ered. For each included species i this results in a balance equa- 
tion of the form: 



(ill 1 v~i v~~i 

-37 = Zj2j n i nkR j ki + Zj 

' i 1 

X& + XX 



j k 
( 

-H; 



V / 



I J 



(4) 



Here n, denotes the density of species i. The first two terms cover 
all formation processes while the last two terms account for all 
destruction reactions. Rm is the reaction rate coefficient for the 
reaction Xj + X* — > X,- + ... (X stands for species X), £7 is 
the local photo-destruction rate coefficient for ionization or dis- 
sociation of species X,- + hv — » X; + either by FUV photons 
or by cosmic ray (CR) induced photons, and £/,• is the local for- 
mation rate coefficient for formation of X, by photo-destruction 
of species X;. For a stationary solution one assumes drij/dt - 0, 
while non-stationary models solve the differential equation 
in time. The chemical network is a highly non-linear system of 
equations. Hence it is not self-evident that a unique solution ex- 
ists a t all, multiple solutio n may b e possible as demonstrat ed e.g. 
bv lLe Bourlot et alj d 19931) and lBoger & Sternberg! d2006l) . 

They showed that bistability may occur in gas-phase models 
(neglecting dust chemistry) of interstellar dark clouds in a nar- 
row parameter range of approximately 10 3 crrT 3 > n/f-17 > 
10 2 cirT 3 with the cosmic-ray ionization rate of molecular hy- 
drogen £cr = 10~ 17 £-n s _1 . Within this range the model re- 
sults may depend very sensitively on variations of input pa- 
rameters such as £cr or the H3 dissociative recombination rate. 
To demonstrate this we show the influence of varying ioniza- 
tion rates in Fig. [4] The left panel gives abundance profiles for 
benchmark model Fl (n=10 3 cirT 3 , x - 10) the right panel 
shows a similar model but with higher density (n=10 4 cm -3 ). 
The higher density was chosen to make sure that we are outside 
the bistability regime. The solid lines in both panels are for a 
cosmic ray helium ionization rate of ^cR(He) = 2.5 x 10~ 17 s _1 , 
the dashed lines denote an ionization rate increased by a fac- 
tor four. Different colors denote different chemical species. The 
most prominent differences are highlighted with colored arrows. 
The factor four in ^cr(Hc) results in differences in density up 
to three orders of magnitude in the lower density case! A de- 
tailed analysis shows that the strong abundance transitions occur 
for ^cR(He) > 8 x 10~ 17 s _1 . This highly non-linear behavior 
disappears if we leave the critical parameter range as demon- 
strates in the right panel of Fig. [4] Boger & Sternberg (2006) 
emphasize that this effect is a property of the gas phase chemical 
network, and is damped if gas-grain processes such as grain as- 
sisted recombination of the at omic ions are introduced (see also 
Shalabiea & Green bergT, 1 1995b . They conclude that the bistabil- 
ity phenomenon pr obably does not occ ur in realistic dusty inter- 
stellar clouds while iLe Bourlol d2006h argues that what matters 
for bistability is not the number of neutralisation channels but the 
degree of ionisation and that bistability may occur in interstellar 



clouds. They suggest this could be one of the p ossible reasons o f 
the non detection of 2 by the ODIN satellite dViti et aUl2001l) . 
Yet, another possible explanation for the absence of O2 is freeze- 
out onto dust. However it is clear that bistability is a real prop- 
erty of interstellar gas-phase networks and not just a numerical 
artifact. Furthermore it is important to emphasize that standard 
PDR models may react very sensitively on the variation of input 
parameters (e.g. £cr. the H2 formation rate, the PAH content of 
the model cloud, etc.) and one has to be careful in the interpre- 
tation of surprising model signatures. 

The numerical stability and the speed of convergence may 
vary significantly over different chemical networks. Three major 
questions have to be addressed: 



1 . which species i are to be included? 

2. which reactions are to be considered? 

3. which reaction rate coefficients are to be applied? 

A general answer to question 1 cannot be given, since this de- 
pends on the field of application. In steady state one has to solve 
a system of M nonlinear equations, where M is the number of 
included species, thus the complexity of the problem scales with 
the number of species (oc N 2 ...N 3 ) rather than with the number 
of chemical reactions. Nowadays CPU time is not a major driver 
for the design of chemical networks. Nevertheless, in some cases 
a small network can give similar results as a big network. Several 
studies have shown that very large networks may include a 
surprisingly large number of 'unimportant' reactions, i.e. reac- 
tions that may be removed from the network without changing 
the chemical structure significantly ( Markwick -Kempeil [2005; 
Wakelam et all l2005al) . It is more important to identify cru- 
cial species not to be omitted, i.e. species that dominate the 
chemical structure under certain conditions. A well known ex- 
ample is the importance of sulfur for the formation of atomic 
carbon at intermediate Ay where the charge transfer reaction 
S + C + — » C + S + constitutes an additional production chan- 
nel f or atomic carbon, visible in a second rise in the abundance 
of C (Sternberg & Dalgarno, 1995). In these benchmarking cal- 
culations, sulfur was not included in order to minimize model 
complexity, in spite of its importance for the PDR structure. 

Regarding question 2 a secure brute force approach would be 
the inclusion of all known reactions involving all chosen species, 
under the questionable assumption that we actually know all im- 
portant reactions and their rate coefficients. This assumption is 
obviously invalid for grain surface reactions and gas-grain in- 
teractions such as freeze-out and desorption. It is important not 
to create artificial bottlenecks in the reaction scheme by omit- 
ting important channels. The choice of reaction rate coefficients 
depends on factors like availability, accuracy, etc.. A number 
of comprehensive datab ases of rate co e fficients is av ailable to- 
day, e.g. NSM/OHIO dWakelam et al.L 12001 |20Q5bh. UM IST 
(Mill anFarquhar. & Willacvl Il997t ILe Teuff et all |2000), and 



Meudon dLe Bourlot et aU 1 19931) . which collect the results from 



many different references, both theoretical and experimental. 

An example for the importance of explicitly agreeing on the 
details of the computation of the reaction rate is the reaction: 



C + H 2 -> CH + H 



(5) 



It has an activation e n ergy barrier of 11700 K 
(Millar. Farquh ar7& Willacvl 1 1997b . effectively reducing 
the production of CH molecules. If we include vibrationally 
excited W 2 into the chemical network and assume that reaction 
(O has no activation energy barrier for reactions with Hi we 
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Fig. 4. The influence of the cosmic ray ionization rate on the chemical structure of a model cloud. The left panel shows results for 
Model Fl (n=10 3 cm 3 , x = 10), the right panel gives results for 10 times higher densities (n=10 4 crrT 3 , x - 10). The solid lines 
give the results for a cosmic ray ionization rate of Helium, enhanced by a factor 4, the dashed lines are for the lower ionization rate. 
The different colors denote different chemical species. The most prominent differences are highlighted with colored arrows. 
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Fig. 5. Comparison between model codes with (dashed line) and 
without (solid line) excited molecular hydrogen, W T The abun- 
dance profile of CH is plotted for both models against Ay.eff. 
Benchmark model F3 has a high density (n = 10 5 5 cirT 3 ) and 
low FUV intensity (x = 10). 



obtain a significantly higher production rate of CH as shown in 
Figure [5]Even this approach is a rather crude assumption, but it 



demonstrates the importance of explicitly agreeing on how to 
handle the chemical calculations in model comparisons. 

Another example is the formation of C in the dark cloud part 
of a PDR, i.e. at values of Ay > 5. A possible formation chan- 
nel for atomic carbon is the dissociat ion of CO by secon dary 
UV photons, induced by cosmic rays dLe Teuff et all |2000|) . In 
the outer parts of the PDR the impinging FUV field dominates 
the dissociation of CO, but for high Ay the FUV field is effec- 
tively shielded and CR induced UV photons become important. 
For CO, this process depends on the l evel population of C O, and 
therefore is temperature dependent dGredel et all 1 1987b . how- 
ever this temp erature depe n dence is often ignored. The reaction 
rate given by Gred el etaT] 0987) has to be corrected by a fac- 
tor of (r/300/T) 1 ' 17 effectivel y reducing the dissoc iation rate for 
temperatures below 300 K dLe Teuff et all 120001) . In Figure |6] 
we plot the density profile of atomic carbon for an isothermal 
benchmark model with temperature T = 50 K. The solid line 
represents the model result for an uncorrected photo-rate using 
the average reaction rate for T = 300 K, compared to the results 
using the rate corrected for T=50 K by (50/300) 117 , given by the 
dashed curve. 



3.1 .3. Heating and Cooling 

To determine the local temperature in a cloud, the equilibrium 
between heating and cooling has to be calculated. The heating 
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Fig. 6. The density profile of atomic carbon for the benchmark 
model F2 (low density, high FUV, T=const=50K, as discussed 
in §|H). The solid curve results from a constant dissociation by 
CR induced secondary photons (implicitely assuming T=300K), 
the dashed curve shows the influence of a temperature depen- 
dent dissociation, i.e. the corresponding dissociation rate was 
corrected by a factor of (77300#) 117 with T=50K. 



rates mainly depend on the H2 formation rate, the electron 
density, the grain size distribution, grain composition, and H2 
treatment (i.e. two-line approximation vs. full ro-vib model), 
while the cooling rates are dominantly influenced by the 
abundance of the main cooling species and the dust opacity 
in the FIR. Table [TJ gives an overview of the most important 
heating and cooling processes. Most of them can be modelled at 
different levels of detail. This choice may have a major impact 
on the model results. One example is the influence of PAHs 
on the efficiency of the photoelectric heating, which results 
in a significantly higher temperature e.g. at t he surface of the 
mode l cloud if PAHs are taken into account. Bake s & Tielensl 
(1994) give convenient fitting formulas for the photoelectric 
heating. Another important case is the collisional de-excitation 
of vibrationally excited H2. A detailed calculation of the level 
population shows that for temperatures above 800 K the lower 
transitions switch from heating to cooling. This imposes a 
significant influence on the net heating from H2 vibrational 
de-excitation. When using an approximation for the heating rate 
it is i mportant to account for this cooling effect dRollig et al.L 
2006). The cooling of the gas by line emission depends on 



Table 1. Overview over the major heating and cooling processes 
in PDR physics 



heating 


cooling 


photoelectric heating (dust & PAH) 


[CII] 158yL/m 


collisional de-excitation of vib. excited H2 


[OI] 63, 145^ 


H2 dissociation 


[CI] 370, 610 yum 


H 2 formation 


[Sill] 35 yum 


CR ionization 


COJFO, OH, H, 


gas-grain collisions 


Ly a, [OI], [Fell] 


dissipation of turbulence 


gas-grain collisions 



the atomic and molecular constants as well as on the radiative 
transfer. A common approximation to the radiative transfer 
probl em is by assuming escape probabili ties for the coolin g 
lines dde Jong et all 1 19801: IStutzkj 1 1984k IStorzer et all Q996). 
The excitation temperature at any point can be computed by 
balancing the collisional excitation and the photon escape prob- 
ability. The local escape probability is obtained by integrating 
exp(-T v ) over An. In the escape probability approximation it 
is now assumed that the radiative interaction region is small 
enough so that the optical depth in each direction is produced 
by molecules with the same excitation temperature. Then the 
excitation problem becomes a local one. The [OI] 63/im line 
may also become very optically thick and can act both as heating 
and cooling contribution. Under certain benchmark conditions 
(low density, constant temperature T gas - 50 K) the [OI] 
63/im line even showed weak maser behavior ( see data plots 
at |http : / /www . phi . uni-koeln . de/pdr-comparison| l. 
Collisions between the gas particles and the dust grains also 
contribute to the total heating or cooling. 

3.1.4. Grain Properties 

Many aspects of PDR physics and chemistry are connected to 
dust properties. We will give only a short overview of the impor- 
tance of dust grains in the modeling of PDRs. Dust acts on the 
energy balance of the ISM by means of photoelectric heating; it 
influences the radiative transfer by absorption and scattering of 
photons, and it acts on the chemistry of the cloud via grain sur- 
face reactions, e.g. the formation of molecular hydrogen and the 
depletion of other species. One distinguishes three dust compo- 
nents: PAHs, very small grains (VSGs) and big grains (BGs). 

The pro perties of big grains have been reviewed recently by 
lDraind(l2003L and references therein). The first indirect evidence 
for the presence of VSGs in the ISM was presented by Andriesse 
(I1978I) in the case of the M17 PDR. The dust grains themselves 
consist of amorphous silicates and carbonaceous material and 
may be covered with ice mantles in the denser and colder parts 
of the ISM. For details of the composition o f grains and their ex- 
tinction due to scattering and absorption see lLi & Draind 12002) 
and references therein. 

The influence and proper treatment of electron densities to- 
gether with grain ionization and recombination is still to be an- 
alyzed. Not only the charge of dust and PAHs but also the scat- 
tering properties are still in discussion ( Weing artner & DraineL 
2001). This heavily influences the model output, e.g. the inclu- 
sion of back-scattering significantly increases the total H2 photo- 
dissociation rate at the surface of the model cloud compared to 
calculations with pure forward scattering. 

3.1 .5. Radiative Transfer 

The radiative transfer (RT) can be split into two distinct wave- 
length regimes: FUV and IR/FIR. These may also be labeled 
as 'input' and 'output'. FUV radiation due to ambient UV field 
and/or young massive stars in the neighborhood impinges on the 
PDR. The FUV photons are absorbed on their way deeper into 
the cloud, giving rise to the well known stratified chemical struc- 
ture of PDRs. In general, reemission processes can be neglected 
in the FUV, considerably simplifying the radiative transfer prob- 
lem. Traveling in only one direction, from the edge to the inside, 
the local mean FUV intensity can usually be calculated in a few 
iteration steps. In contrast to the FUV, the local FIR intensity is 
a function of the temperature and level populations at all posi- 
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tions due to absorption and reemission of FIR photons. Thus a 
computation needs to iterate over all spatial grid points. A com- 
mon simplifying approximation is the spatial decoupling via the 
escape probability approximation. This allows to substitute the 
intensity dependence by a dependence on the relevant optical 
depths, ignoring the spatial variation of the source function. The 
calculation of emission line cooling then becomes primarily a 
problem of calculating the local excitation state of the particular 
cooling specie s. An overview of NLTE ra diative transfer meth- 
ods is given bv lvan Zadelhoffetindljoi 



4. Description of the Benchmark Models 

4.1. PDR Code Characteristics 

A total number of 11 model codes participated in the PDR 
model comparison study during and after the workshop in 
Leiden. Table|2]gives an overview of these codes. The codes are 
different in many aspects: 



- finite and semi-infinite plane-parallel and spherical geome- 
try, disk geometry 

- chemistry: steady state vs. time-dependent, different chemi- 
cal reaction rates, chemical network 

- IR and FUV radiative transfer (effective or explicitly wave- 
length dependent), self- and mutual shielding, atomic and 
molecular rate coefficients 

- treatment of dust and PAHs 

- treatment of gas heating and cooling 

- range of input parameters 

- model output 

- numerical treatment, gridding, etc. 

This manifold in physical, chemical and technical differences 
makes it difficult to directly compare results from the differ- 
ent codes. Thus we tried to standardize the computation of the 
benchmark model clouds as much as possible. This required all 
codes to reduce their complexity and sophistication, often be- 
yond what their authors considered to be acceptable, consid- 
ering the actual knowledge of some of the physical processes. 
However as the main goal of this study was to understand why 
and how these codes differ these simplifications are acceptable. 
Our aim was not to provide a realistic model of real astronomical 
objects. The individual strengths (and weaknesses) of each PDR 
code are briefly summarized in the Appendix and on the website: 
|http : / /www. phi . uni-koeln . de/pdr-comparison| . 

4.2. Benchmark Frame and Input Values 

A total of 8 different model clouds were used for the benchmark 
comparison. The density and FUV parameter space is covered 
by accounting for low and high densities and FUV fields under 
isothermal conditions, giving 4 different model clouds. In one set 
of models the complexity of the model calculations was reduced 
by setting the gas and dust temperatures to a given constant value 
(models F1-F4, 'F' denoting a fixed temperature), making the 
results independent of the solution of the local energy balance. 
In a second benchmark set, the thermal balance has been solved 
explicitly thus determining the temperature profile of the cloud 
(models V1-V4, 'V denoting variable temperatures). Table [3] 
gives an overview of the cloud parameter of all eight benchmark 
clouds. 



Table 3. Specification of the model clouds that were computed 
during the benchmark. The models F1-F4 use constant gas and 
dust temperatures, while VI -V4 have their temperatures calcu- 
lated self consistently. 





Fl 






F2 






T=50K 






T=50K 




n = 


10 3 cm" 3 ,^ = 


10 


n = 


10 3 cm- 3 ,^ = 


10 5 




F3 






F4 






T=50K 






T=50K 




n = 


10" m 3 ,^ = 


10 


n = 


10" crrT 3 ,^ = 


10 5 




VI 






V2 






T=variable 






T=variable 




n = 


10 3 crrr 3 ,,y = 


10 


n = 


10 3 crrr 3 ,^ = 


10 5 




V3 






V4 






T=variable 






T=variable 




n = 


10" cmr\ X = 


10 


n = 


10" crrT 3 ,^ = 


10 5 



4.2.1. Benchmark Chemistry 

One of the crucial steps in arriving at a useful code comparison 
was to agree on the use of a standardized set of chemical species 
and reactions to be accounted for. For the benchmark models we 
only included the four most abundant elements H, He, O, and C. 
Additionally only the species given in Tab. [4] are included in the 
chemical network calculations: 



Table 4. Chemical content of the benchmark calculations. 



Chemical species in the models 

H, H + , H2, H9, H^ 

O, + , OFF, OH, O,, O}, H 2 0, H 2 + , H 3 + 
C, C + , CH, CH+, CH 2 , CH+, CH 3 , 
CH3, CH 4 , CH+, CH+, CO, CO + ,HCO + 
He, He + , e" 



The chemical reactio n rates ar e taken from the 
UMIST99 database dLe Teuff et al.L l2000h together 
with some corrections suggested by A. Sternberg. 
The complete reaction rate file is available online 
( |http : / /www . phi . uni-koeln . de/pdr- comparison} . 
To reduce the overall modeling complexity, PAHs were ne- 
glected in the chemical network and were only considered 
for the photoe lectric heating (photo electric heating efficiency 
as given by iBakes & TielensL 1 19941) in models V1-V4. Codes 
which calculate time-dependent chemistry used a suitably long 
time-scale in order to reach steady state (e.g. UCL_PDR used 
100 Myr). 

4.2.2. Benchmark Geometry 

All model clouds are plane-parallel, semi-infinite clouds of con- 
stant total hydrogen density n = n(H) + 2 n(H2). Spherical codes 
approximated this by assuming a very large radius for the cloud. 

4.2.3. Physical Specifications 

As many model parameters as possible were agreed upon at the 
start of the benchmark calculations, to avoid confusion in com- 
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Model Name 


Authors 


Cloudy 


G. J. Ferland, P. van Hoof, N. P. Abel, G. Shaw (Ferland et al., 1998; Abel et al., 2005; Shaw et al., 2005) 


COSTAR 


I. Kamp, F. Bertoldi, G.-J. van Zadelhoff (Kamp & Bertoldi, 2000; Kamp & van Zadelhoff, 2001) 


HTBKW 


D. Hollenbach, A.G.G.M. Tielens, M.G. Burton, M.J. Kaufman, M.G. Wolfire 
( lielens & Hollenbach, 19oj; Kaurman et al., 1999; Wolnre et al., zl)U3) 


KOSMA-T 


H. Storzer, J. Stutzki, A. Sternberg (Storzer et al., 1996), B. Koster, M. Zielinsky, U. Leuenhagen 
rsenscn et ai. (zuujj.Koiiig et ai. (zuuoj 


Lee96mod 


H.-H. Lee, E. Herbst, G. Pineau des Forets, E. Roueff, J. Le Bourlot, O. Morata (Lee et al., 1996) 


Leiden 


J. Black, E. van Dishoeck, D. Jansen and B. Jonkheid 

(Black & van Dishoeck, 1987; van Dishoeck & Black, 1988; Jansen et al., 1995) 


Mei jerink 


R. Meijerink, M. Spaans (Meiierink & Spaans, 2005) 


Meudon 


J. Le Bourlot, E. Roueff, F. Le Petit (Le Petit et al., 2005, 2002; Le Bourlot et al., 1993) 


Sternberg 


A. Sternberg, A. Dalgarno (Sternberg & Dalgarno, 1989, 1995; Boger & Sternberg, 2005) 


UCL.PDR 


S. Viti, W.-F. Thi, T. Bell (Tavlor et al., 1993; Papadopoulos et al., 2002; Bell et al., 2005) 



paring model results. To this end we set the most crucial model 
parameters to the following values: the value for the standard 
UV field was taken as x - 10 an d 10 5 times the Draine (1978) 
field. For a semi-infinite plane parallel cloud the CO dissociation 
rate at the cloud surface for^- = 10 should equal 1CL 9 s _1 , using 
that for optically thin conditions (for which a point is exposed to 
the full An steradians, as opposed to 2n at the cloud surface) the 
CO dissociation rate is 2 x 1CL 10 s~' in a unit Draine field. The 
cosmic ray H ionization rate is assumed to be ( — 5 x 1CL 17 s" 1 
and the visual extinction Ay = 6.289 x 10 _22 A^e.tot- If the codes 
do not explicitly calculate the H2 photo-dissociation rates (by 
summing over oscillator strengths etc.) we assume that the unat- 
tenuated H2 photo-dissociation rate in a unit Draine field is 
equal to 5.18 x 1CL 11 s , so that at the surface of a semi- 
infinite cloud for 10 times the Draine field the H2 dissociation 
rate is 2.59 x 10~ 10 s" 1 (numerical values from Sternberg. 
See § 15.11 for a discussion on H2 dissociation rates). For the 
dust attenuation factor in the H2 dissociation rate we assumed 
exp(-£Ay) if not treated explicitly wavelength dependent. The 
value k = 3.02 is representative for the effective opacity in the 
912-1120 A range (for a specific value of Ry « 3). We use a 
very simple H2 for mation rate coefficient R ' = 3 x 10" 18 r 1/2 = 
2.121x 17 cm 3 s 1 (ISternberg & Dalearnol [19951) at T = 50 K, 
assuming that every atom that hits a grain sticks and reacts to H2. 
A summary of the most important model parameters is given in 
Table 



5. Results 

In the following section we give a short overview of the up 
to date results of the PDR model comparison. The names of 
the model codes are printed in typewriter font (e.g. COSTAR). 
We will refer to the two stages of the benchmarking results 
by pre- and post-benchmark, denoting the model results 
at the beginning of the comparison and at its end respec- 
tively. All pre- and post-benchmark results are posted at 
|http : //www . phi . uni-koeln . de/pdr-comparison| One 
model from the initial 12 participating model was left out in 
the post-benchmark plots because the authors cou ld not attend 
the w orkshop. In addition, the KOSMA-t models (Rollig Tt all 
2006) and the models by Bensch, which participated in the 
comparison as seperate codes, have been merged to a single 
set (labeled KOSMA-t) as they are variants on of the same 
basic model which do not differ for the given benchmarking 



parameter set, and consequently give identical results. To 
demonstrate the impact of the benchmark effort on the results of 
the participating PDR codes we plot the well known CI C + / CO 
transition for a typical PDR environment before and after the 
changes identified as necessary during the benchmark in Fig. [7] 
The photo-dissociation of carbon mon oxide is thought to be 
well understood for almost 20 years (Ivan Dishoeck & B lack. 
1988). Nevertheless we see a significant scatter for the densities 
of C + , C, and CO in the top plot of Fig. [7] The scatter in the 
pre-benchmark rates is significant. Most deviations could be 
assigned to either bugs in the pre-benchmark codes, misunder- 
standings, or to incorrect geometrical factors (e.g. 2n vs. 4n). 
This emphasizes the importance of this comparative study to 
establish a uniform understanding about how to calculate even 
these basic figures. Despite the considerable current interest 
because of, e.g. SPITZER results, we do not give the post- 
benchmark predictions for the H2 mid-IR and near IR lines (or 
the corresponding Boltzmann diagram). Only a small fraction 
of the participating codes is able to compute the detailed H2 
population and emission, and the focus of this analysis is the 
comparison between the benchmark codes. 



Table 5. Overview of the most important model parameter. All 
abundances are given w.r.t. total H abundance. 



Model Parameters 



A He 


0.1 


elemental He abundance 


A 


3 x 10~ 4 


elemental O abundance 


A c 


1 x 10~ 4 


elemental C abundance 


fcs 


5 x 10~ 17 s- [ 


CR ionization rate 


A v 


6.289 x W- 12 N Hloul 


visual extinction 


T UV 


3.02A V 


FUV dust attenuation 


n 


1 km s~' 


Doppler width 




5 x 10- 18 • $ s-' 

3 x 10- 18 r'^ cm 3 s-' 


H2 dissociation rate 


R 


H 2 formation rate 




50 K 


gas temperature (for F1-F4) 




20 K 


dust temperature (for F1-F4) 


n 


10 3 , 10" cm" 3 


total density 


X 


10, 10 5 


FUV intensity w.r.t. 



Draine (1978) field 
(i.ejf = 1.71 G ) 
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Fig. 7. Model Fl (n=10 3 cm 3 , x — 10): Comparison between the density profiles of C + (top), C (middle), and CO (bottom) before 
(top) and after (bottom) the comparison study. The vertical lines indicate the code dependent scatter. For C and CO they indicate the 
depths at which the maximum density is reached, while for C + they indicate the depths at which the density dropped by a factor of 
10. Dashed lines indicate pre-benchmark results, while solid lines are post-benchmark. 
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V.eff v.eff V.eff 

Fig. 8. Model Fl (n=10 3 cm -3 , \ — 10): The photo-dissociation rates of H2 ( left column), of CO (middle column) and the photo- 
ionization rate of C ( right column) after the comparison study. 



5. 1. Models with Constant Temperature F1-F4 

The benchmark models Fl to F4 were calculated for a fixed 
gas temperature of 50 K. Thus, neglecting any numerical issues, 
all differences in the chemical structure of the cloud are due 
to the different photo-rates, or non-standard chemistry. Some 
PDR codes used slightly different chemical networks. The code 
Sternberg uses the standard chemistry with the addition of 
vibrational excited hydrogen and a smaller H-H2 formation net- 
work. The results by Cloudy were obtained with two different 
chemical setup s: The pre-benchmark c hemi stry had the chem- 
ical network of iTielens & Hollenbachl d!985l) . The post bench- 
mark results use the corrected UMIST database. Cloudy also 
used a different set of radiative recombination coefficients for 
the pre-benchmark cal culations which w ere the major source for 
their different results dAbel et all I2005I) . The carbon photoion- 
ization and radiative recombination rates are very sensitive to ra- 
diative transfer and hence to dust properties. The dust properties 
in Cloudy are different from what is implicitly assumed in the 
UMIST fits. Cloudy' s post-benchmark results are achieved 
after switching to the benchmark specifications. After the switch 
they agree very well with the other codes. In Fig. |7]we present 
the pre- and post-benchmark results for the main carbon bearing 
species C + , C, and CO. To emphasize the pre-to-post changes 
we added several vertical marker lines to the plots. For C and 
CO they indicate the depths at which the maximum density is 
reached, while for C + they indicate the depths at which the 
density has dropped by a factor of 10. Dashed lines indicate 
pre-benchmark results, while solid lines are post-benchmark. In 
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Fig. 9. Model Fl (n=10 3 cm" 3 , x = 10) The H " H 2 transition 
zone after the comparison study. Plotted is the number density 
of atomic and molecular hydrogen as a function of A v , e tf- The 
vertical lines denote the range of the predicted transition depths 
for pre- and post-benchmark results (dashed and solid lines re- 
spectively). 



the pre-benchmark results the code dependent scatter for these 
depths is A Av, e ff ~ 2 - 4 and drops to A Ay.eff ~ 1 in the post- 
benchmark results. 
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In the post-benchmark results, the Leiden and UCL_PDR 
models show a slightly different behavior. The predicted peak 
depth of C is somewhat smaller than for the other codes. The 
peak C density of UCL_PDR is roughly 50% higher than in the 
other codes. A comparison with the photo-ionization of C shown 
in Fig.[8]confirms that a slightly stronger shielding for the ioniza- 
tion of C is the reason for the different behavior of C and C + . The 
dark cloud densities for C + , C, and CO agree very well, except 
for a somewhat smaller C + density in the Lee 9 6mod results. 

In Fig. [8] we plot the post-benchmark photo-rates for disso- 
ciation of H2 (left column) and CO (middle column) and for the 
ionization of C (right column), computed for model Fl . Even for 
this simple model there are some significant differences between 
the models in the various rates. In the pre-benchmark results, 
several codes calculated different photo-rates at the edge of the 
model cloud, i.e. for very low values of Ay, e fr- Some codes cal- 
culated surface photo-dissociation rates between 4-5 x 1CT 10 s _1 
compared to the expected value of 2.59x 10~ 10 s _1 . Most of these 
deviations were due to exposure to the full An steradians FUV 
field instead the correct 2n, but also due to different effects, like 
the FUV photon back-scattering in the Meudon results. The pre- 
benchmark rates of KOSMA-r were shifted toward slightly lower 
values of Ay because of an incorrect scaling between Ay and 
Ay >e ff and an incorrect calculation of the angular averaged photo- 
rate (the model features a spherical geometry with isotropic FUV 
illumination). The post-benchmark results (Fig. [8]) show that 
most deviations have been corrected. The remaining offset for 
the Meudon result is due to the consideration of backscattered 
FUV photons, increasing the local mean FUV intensity. The pre- 
to post-benchmark changes for the photo-rates of CO and C are 
even more convincing (see online archive). The post-benchmark 
results are in very good agreement except for some minor differ- 
ence, e.g. UCL_PDR' s photo-ionization rate of C showing some 
deviation from the main field. 

The depth-dependence of the H2 photo-dissociation rate is 
reflected in the structure of the H-H2 transition zone. Fig. [9] 
shows the densities of atomic and molecular hydrogen after the 
benchmark. The vertical lines denote the minimum and max- 
imum transition depths before (dashed) and after the bench- 
mark (solid). In the pre-benchmark results the predicted tran- 
sition depth ranges from 0.08 Ay e ff to 0.29 Ay.eff- In the post- 
benchmark results this scatter is reduced by more than a factor 
of 3. Sternberg gives a slightly smaller H density in the dark 
cloud part. In this code, cosmic ray (CR) destruction and grain 
surface formation are the only reactions considered in the calcu- 
lation of the H2 density. The other codes use additional reactions. 
The reactions: 

Hj + H 2 — > + H (/fc = 2.08x 10 -9 cm 3 s -1 ) 
H 2 + CH+ -> CH3 + H (Jt = 1.6 X 10" 9 cm 3 s _1 ) 

contribute to the total H density at high Ay, e ff- This results in 
a somewhat higher H density as shown in Fig. [9] The Meudon 
model gives a slightly smaller H2 density at the edge of the cloud 
than the other codes. This is due to the already mentioned higher 
photo-dissociation rate of molecular hydrogen applied in their 
calculations. 

The model Fl may represent a typical translucent cloud 
PDR, e.g., the line of sight toward HD 147889 in Ophiuchus 
dLiseau et al. LI1999I) . The low density and FUV intensity condi- 
tions emphasize some effects that would be hard to notice oth- 
erwise. This includes purely numerical issues like gridding and 
interpolation/extrapolation of shielding rates. These differences 
explain why the various codes still show some post-benchmark 



scatter. We relate differences in the predicted abundances to the 
corresponding rates for ionization and dissociation. 

Since most of the codes use the same chemical network and 
apply the same temperature, the major source for remaining de- 
viations should be related to the FUV radiative transfer. To study 
this we present some results of benchmark model F4 featuring 



a density n = 10 



and a FUV intensity x - 10 



order to enhance any RT related differences and discuss them 
in more detail. Fig. [10] shows the post-benchmark photo-rates 
for the model F4. The higher unshielded H2 photo-rate in the 
Meudon results, already visible in model Fl (Fig. [8]l is now sig- 
nificantly enhanced due to the increased FUV flux. Meudon, 
as well as Cloudy, Leiden and Sternberg, treat the hy- 
drogen molecule by calculating the local level population and 
determining the photo-dissociation rate by integrating each ab- 
sorption line over the absorption cross section and summing 
over all absorption lines. Meudon, Cloudy, and Leiden in- 
tegrate the line profile over the attenuated spectrum, in order 
to account for line overlap effects, while Sternberg treats 
each line seperately, neglecting line overlap. Most other codes 
just assume that the photodissociation scales with the incident 
radiation field, neglecting any influence from varying H2 level 
populations. One reason for the different H2 photo-rate is a dif- 
ferent local mean FUV intensity, caused by backscattered pho- 
tons. However, this should only account for approximately 10% 
of the increased dissociation rate. The remaining differences are 
due to different treatment of H2. Either they use different equa- 
tions, e.g. full ro-vib resolution in Meudon and Sternberg vs. 
only vib. population in KOSMA-r, or they apply different molec- 
ular d a ta. Sternberg uses data f rom ISternberg & Dalgarnol 
(119891) : ISternberg & Neufeldl d!999l). Meudon uses co llisional 
data from Flower (ll997lll998l) : lFlower & Rouefj(fl9 99) and as- 
sociated papers, and radiative data from lAbgrall et al.l(l2000l) . in- 
cluding dissociation efficiencies. These different data sets result 
in: 

1. Excited rotational states are much more populated in 
Meudon' s results than in Sternberg 

2. Dissociation from an excited rotational level increases much 
faster with J in Meudon' s data. 

Both effects lead to dissociation probabilities that differ by 2-3 in 
case of Model F4. Due to the structure of the code these features 
could not be turned off in Meudon results. 

The photo-rates for CO and C are in very good accord, but 
we notice a considerable spread in the shielding behavior of the 
hydrogen photo-rate. This spread is due to the particular imple- 
mentation of H2 shielding native to every code, by either using 
tabulated shielding functions or explicitly calculating the total 
cross section at each wavelength. The different photo-rates di- 
rectly cause a different H-H2 transition profile, shown in the top 
panel of Fig. QT| The low molecular hydrogen densities in the 
Meudon and Cloudy models are again due to the higher H2 
photo-dissociation rate. Sternberg' s slightly lower H2 abun- 
dance at the edge of the cloud is consistent with the marginally 
higher, unshielded H2 photo-dissociation rate, seen in the top 
plot in Fig. [TO] The Mei jerink code shows the earliest drop 
in the photo-rate, reflected by the corresponding increase in the 
H2 density. The qualitatively different H2 profile in KOSMA-t 
is most likely due to the spherical geometry in the code. Again 
Sternberg produces slightly smaller H densities for high val- 
ues of Ayeff. Since Sternberg does not consider additional 
reactions for the H/H2 balance its H density profile is the only 
one not showing the slight kink at Ay, e ff ~ 2. ..3. These devia- 
tions do not significantly change the total column density of hy- 
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Fig. 10. Model F4 (n=10 5 5 cm 3 , x — 10 5 ): The post-benchmark photo-dissociation rates of H2 (left column), of CO (middle 
column) and the photo-ionization rate of C (right column) (upper plot). 



drogen. Hence the impact on any comparison with observational 
findings is small. Nevertheless one would expect that under the 
standardized benchmark conditions all codes produce very sim- 
ilar results, yet we note a considerable spread in hydrogen abun- 
dances for Ay.eff > 2. This again emphasizes how complex and 
difficult the numerical modeling of PDRs is. The bottom panel 
in Fig. [TTJ shows the density profiles of C + , C, and CO. Here, 
the different codes are in good agreement. The largest spread 
is visible for the C density between Av, e ff ~ 3. ..6. The results 
for C + and CO differ less. Lee96mod' s results for C + and C 
show a small offset for Av, e ff > 6. They produce slightly higher 
C abundances and lower C + abundances in the dark cloud part. 
The different codes agree very well in the predicted depth where 
most carbon is locked up in CO (Ay, e ff ~ 3. 5. ..4. 5). This range 
improved considerably compared to the pre-benchmark predic- 
tions Of Ay.eff ~ 3. ..8. 

The results from models F1-F4 clearly demonstrate the 
importance of the PDR code benchmarking effort. The pre- 
benchmark results show a significant code-dependent scatter. 
Although many of these deviations have been removed during 
the benchmark activity, we did not achieve identical results with 
different codes. Many uncertainties remained even in the isother- 
mal case, raising the need for a deeper follow up study. 

5.2. Models with Variable Temperature V1-V4 

In the benchmark models V 1 - V4 the various codes were required 
to also solve the energy balance equations in order to derive the 
temperature structure of the model clouds. This of course intro- 
duces an additional source of variation between the codes. The 



chemical rate equations strongly depend on the local tempera- 
ture, hence we expect a strong correlation between temperature 
differences and different chemical profiles of the model codes. 
As a consequence of a differing density profile of e.g. CO and 
H2 we also expect different shielding signatures. We will restrict 
ourselves to just a few exemplary non-isothermal results because 
a full analysis of the important non-isothermal models goes be- 
yond the scope of this paper. To demonstrate the influence of a 
strong FUV irradiation we show results for the benchmark model 
V2 with n = 10 3 cirT 3 , and^ = 10 5 in Figs. ll2H16l The detailed 
treatment of the various heating and cooling processes differs 
significantly from code to code. The only initial benchmark re- 
q uirements was to treat th e photoelectric (PE) heating according 
to lBakes & Tielensl(ll994h . On one hand, this turned out to be not 
strict enough to achieve a sufficient agreement for the gas tem- 
peratures, on the other hand it was already too strict to be eas- 
ily implemented for some codes, like Cloudy, which calculates 
the PE heating self-consistently from a given dust composition. 
This demonstrates that there are limits to the degree of standard- 
ization. The calculation of the dust temperature was not stan- 
dardized and varies from code to code. Since Lee96mod only 
accounts for constant temperatures, their model is not shown in 
the following plots. We only plot the final, post-benchmark sta- 
tus. 

In Fig. [12] we show the gas temperature over Av, e ff- The gen- 
eral temperature profile is reproduced by all codes. Even so we 
note considerable differences between different codes. The de- 
rived temperatures at the surface vary between 1600 and 2500 K. 
For low values of Ay. e ff the heating is dominated by PE heating 
due to the high FUV irradiation, and the main cooling is pro- 
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Fig. 11. Model F4 (n=10 cm 3 , x = 10 5 ): The upper panel shows the post-benchmark results for the H and H2 densities. The lower 
panel shows the post-benchmark density profiles of C + , C, and CO. The vertical gray lines in both panels indicate the pre-to-post 
changes. 
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Fig. 12. Model V2 (n=10 3 cirT 3 , x = 10 5 ): Th e plot shows the 
post-benchmark results for the gas temperature. 



vided by [OI] and [CII] emission. It is interesting, that the dom- 
inant cooling line is the [OI] 63/mi line (cf. Fig. [T6j left plot), 
although its critical density is two orders of magnitude higher 
than the local density (« cl s5x 10 5 cirT 3 ). The highest surface 
temperature is calculated by Leiden, while Meudon computes 
the lowest temperature. The bulk of models gives surface tem- 
peratures near 1900 K. All models qualitatively reproduce the 
temperature behavior at higher values of Ay.eff and show a min- 
imum temperature of 10 K between Av, e ff ~ 5... 10, followed 
by a subsequent rise in temperature. The only relevant heating 
contribution at Ay.eff > 5 comes from cosmic ray heating, which 
hardly depends on Av, e fl ■ At Av, e ff > 4, the dominant cooling is 
by [CI] fine structure emission. This is a very efficient cooling 
process and the temperature reaches its minimum. At Ay.eff =10 
the atomic carbon density rapidly drops and CO cooling starts 
to exceed the fine structure cooling (cf. abundance profiles in 
Fig- El. However, cooling by CO line emission is much less ef- 
ficient, especially at these low total densities, and thus the tem- 
perature increases again. 

For the bulk of the cloud the heating contribution by H2 
vibrational deexcitation is negligible compared to photoelec- 
tric heating. Only Meijerink and Leiden predict compa- 
rable contributions from both processes. Unfortunately, the ex- 
act treatment of this process was not standardized and depends 
very much on the d etailed implementation (eg. the two-le vel 
approximation fro m iBurton. Hollenbach. & Tielensl dl990h or 
Rolli g et ail d2006l) vs. the solution of the full H2 problem like in 
Meudon, Cloudy, and Sternberg). Generally the heating 
by H2 vibrational deexcitation depends on the local density and 
the local mean FUV intensity, and should thus decrease at large 
values of Av, e ff and dominate the heating for denser clouds. 

At Av, e ff ~ 2.. .3 we note a flattening of the temperature 
curve in many models, followed by a steeper decline some- 
what deeper inside the cloud. This is not the case for HTBKW, 
KOSMA-t, and Sternberg. The reason for this difference is 
the [OI] 63//m cooling, showing a steeper decline for the three 
codes (Fig. [16] left plot). For very large depths, KOSMA-t pro- 
duces slightly higher gas temperatures. This is due to the larger 
dust temperature and the strongest H2 vibrational deexcitation 
heating at Av, e ff > 10. 

In Fig.[l3]we plot the photodissociation rate of H2 (top left), 
the photoioniozation rate of C (top right), and the density of H 
and H2 over Av, e ff (bottom). Meudon' s unshielded dissociation 



rate is by a factor three larger than the median of 2.6 x 10~ 6 s _1 , 
and the Sternberg value of 3.8 x 10 -6 s _1 is slightly larger for 
the same reason as discussed in section 15.11 The depth depen- 
dent shielding shows good agreement between all codes, with 
slight variations. The different model geometry of KOSMA-t 
is reflected in the slightly stronger shielding. Leiden has the 
weakest shielding. Like some of the other codes (see Appendix), 
they account for the detailed H2 problem when calculating the 
photodissociation rate, instead of applying tabulated shielding 
rates. Yet these differences are small, since we are in a parameter 
regime ix/n = 100), where the m ain shielding is dominat ed by 
dust rather than by self shielding (Draine & Bertoldi, 1996). The 
density profiles of H and H2 are in good agreement. The stronger 
photodissociation in Meudon is reflected in their smaller H2 
density at the surface. All other H2 densities correspond well 
to their dissociation rates except for Cloudy, which has a lower 
density at the surface without a corresponding photodissociation 
rate. This is a temperature effect. Cloudy computes relatively 
low surface temperatures which lead to slightly lower H densi- 
ties at the surface. The central densities are also in good accord. 
The different H densities reflect the corresponding temperature 
profiles from Fig. [T2l 

The photoionization rate of C is given in the top right plot 
in Fig. [13] All models are in good agreement at the surface of 
the cloud. Meudon and UCL_PDR drop slightly earlier than the 
bulk of the results. This is also reflected in their C density pro- 
files in Fig. [14] (top right) which incline slightly earlier. Deep 
inside the cloud Sternberg and HTBKW show a steeper de- 
cline compared to the other codes. The agreement for the C + 
profile is also very good. At Av, e ff = 5 the densities drop by a 
factor of 10 and remain constant until they drop at Ay.eff ~ 10. 
This plateau is caused by the increase in C density, compensat- 
ing the FUV attenuating. Leiden' s results show some devia- 
tions for Av >e ff > 10. Their C density remains higher throughout 
the center, causing a slightly different carbon and oxygen chem- 
istry at Av.eff > 10. The calculated O and O2 densities are given 
in Fig. [14] (bottom, right). The dark cloud densities are in very 
good agreement among the models, with some deviations in the 
Leiden values. The O2 profiles show some variations between 
^v.eff ~ 1 and 10 but these are minor deviations especially tak- 
ing the fact that the densities vary over 14 orders of magnitude 
from the outside to the center of the cloud! The differences in 
2 are also reflected in the CO plot (Fig. [HI bottom left). All 
codes produce very similar density profiles and dark cloud val- 
ues. Leiden gives a smaller CO density beyond Ay, e ff = 10. 

In Fig. [15] we plot the total surface brightnesses of the main 
fine-structure cooling lines for the V2 model: [CII] 158 /mi, 
[OI] 63, and 146 /mi, and [CI] 610 and 370 /mi. For the spher- 
ical PDR models, the surface brightness averaged over the pro- 
jected area of the clump is shown. The surface brightness of 
these fine-structure lines is smaller by typically a few 10%, if 
calculated along a pencil-beam toward the clump center as they 
are enhanced in the outer cloud layers. Compared with the pre- 
benchmark results, the spread in Tb has been decreased sig- 
nificantly from almost 3 orders of magnitude to a factor of 3- 
5 for [CII] and [OI]. To explain the differences in Fig. [TBI we 
plot in Fig. Q~6] the radial profiles of the local emissivities of 
[OI] 63/mi and [CI] 310//m. Leiden gives the highest [OI] 
brightnesses and also computes higher local [OI] 63 /mi emis- 
sivities for small values of Av, e ff, shown in Fig. [16] COSTAR, 
with very similar results for the density profile and compa- 
rable gas temperatures, gives much smaller emissivities. The 
reason for these deviations is still unclear. The model depen- 
dent spread in surface brightnesses is largest for the [CI] lines. 



16 



Rollig et al.: A PDR-Code Comparison Study 





Fig. 13. Model V2 (n=10 3 cm 3 , \ - 10 5 ): The post-benchmark photo-dissociation rates of H2 (left column), and the photo- 
ionization rate of C (right column) (upper plot). The lower plots shows the H and H2 densities. 
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Fig. 15. Model V2 (n=10 3 cm 4 , X = 10 5 ): The plot shows the 
post-benchmark surface brightnesses of the main fine-structure 
cooling lines: [CII] 158 fim, [OI] 63, and 146 /mi, and [CI] 610 
and 370 fim. 



HTBKW computes 10 times higher line intensities for the [CI] 
370/im transition than Sternberg. This can be explained as 
follows. Both codes show almost identical column densities and 
abundance profiles of C°, yet the local emissivities are very dif- 
ferent between A v, e ff = 4. ..9 (Fig. |T6l>. Sternberg, together 



with some other codes, compute a local minimum for the cool- 
ing at Av >e ff ~ 6, while the HTBKW, Cloudy, Meijerink, 
and Meudon models peak at the same depth. This can be ex- 
plained as a pure temperature effect, since the codes show- 
ing a [CI] peak compute a significantly higher temperature at 
Av.eff = 6: T(HTBKW)=83 K, T(Sternberg)=10 K. These dif- 
ferent temperatures at the C° abundance peak strongly influ- 
ences the resulting [CI] surface brightnesses. Overall, the model- 
dependent surface temperatures still vary significantly. This is 
due to the additional nonlinearity of the radiative transfer prob- 
lem, which, under certain circumstances, amplifies even small 
abundance/temperature differences. 

5.3. Review of participating codes 

It is not our intent to rate the various PDR model codes. Each 
code was developed with a particular field of application in mind 
and is capable to fulfill its developers expectations. The restric- 
tions artificially posed by the benchmark standards were addi- 
tionally limiting the capacity of the participating model codes. 
Some models encountered for example major numerical diffi- 
culties in reaching a stable temperature solution for the bench- 
mark models V4, mainly caused by the requested H2 formation 
rate of R = 3 x 10 _I8 r^ 2 cm 3 s _1 . This gives reasonable results 
for low temperatures, but diverges for very high temperatures, 
resulting in an unreasonably high H2 formation heating. Other 
codes also show similar numerical problems especially for the 
model V4. This numerical noise vanishes when we apply more 
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Fig. 16. Model V2 (n=10 3 cm 3 , x - 10 5 ): The post-benchmark local emissivities of [OI] 63/im (left column), and [CI] 310/im 
(right column). 



physically reasonable conditions. Nevertheless it was very in- 
structive to study the codes under these extreme conditions. 

Every participating code has its own strengths. The Me u don 
code and Cloudy are certainly the most complex codes in the 
benchmark, accounting for most physical effects by explicit cal- 
culations, starting from the detailed micro-physical processes, 
making the least use of fitting formulae. 

Cloudy was originally develope d to simulate extreme e nvi- 
ronments near accreting black holes dFerland & Reeslfl988h . al- 
though i t has been applied to HII regions, planetary nebulae, and 
the ISM. lFerland et alJ (fT994) describe an early PDR calculation. 
Its capabil ities have been greatly extended over the past sev- 



eral y ears (Ivan Hoof et al 1 12004 lAbel etail 12005b Ishaw et all 
2005). Due to the complexity of the code, it was initially not pos- 
sible to turn off all implemented physical processes as required 
for the benchmark, but during this study they were able to adopt 
all benchmark requirements thus removing all major deviations. 

The codes HTBKW, Leiden, Sternberg and KOSMA-r 
are based on PDR models that began their development 20 years 
ago and have been supported and improved since then. One of 
the main differences between them is the model geometry and il- 
lumination. Plane-parallel geometry and uni-directional illumi- 
nation is assumed in HTBKW, Leiden and Sternberg and 
spherical geometry with an isotropically impinging FUV field 
in KOSMA-r. The chemistry adopted generally in HTBKW is 
the smallest (46 species) compared with Sternberg (78) and 
Leiden/KOSMA-r (variable). Leiden, Sternberg and 
KOSMA-t explicitly solve the H2 problem (full ro-vib level pop- 
ulation) and determine the corresponding shielding by integrat- 
ing all absorption coefficients while HTBKW uses shielding func- 



tions and a single-line approximation for H2. Cloudy is also 
capable of explicitly calculating a fully (v,J) resolved H2 model, 
but this capability was switched off in the final model. Instead 
they used a 3-level approximation there. Leiden and Meudon 
are the only codes in the benchmark explicitly calculating the 
CO shielding, all other codes use shielding factors. HTBKW is ad- 
ditionally accounting for X ray and PAH heating and computes a 
large number of observational line intensities, while Le i den fo- 
cuses on the line emission from the main PDR coolants C + , C, O, 
and CO. However it is possible to couple their PDR output with 
a more sophisticated radiative tran sfer code such as RATRAN 
dHogerheiide & van der Tald , l2000h to calculate emission lines . 
This is also done by KOSMA-t, using O NION dGierens et all 
1 19921) or SimLine dOssenkopf et alll200ll) . COSTAR was devel- 
oped in order to model circumstellar disks, featuring any given 
disk density profile in radial direction and scale height in verti- 
cal direction. It uses uni-directional FUV illumination and can 
treat a surrounding isotropic interstellar FUV field in addition 
to the uni-directional stellar field. It computes a relatively small 
chemical network (48 species) but also accounts for freeze-out 
onto grains and desorption effects. It relies on shielding func- 
tions for H2 and CO and does not calculate observational line 
intensities up to now. Nevertheless most of the COSTAR re- 
sults are in good agreement with the other code results for most 
of the benchmark models, demonstrating that it correctly ac- 
counts for the important PDR physics and chemistry. UCL_PDR 
is a plane-parallel model focused on time-dependent chemistries 
with freeze-out and desorption. Its main features are a fully time- 
dependent treatment - including time-varying density and radi- 
ation profiles - and its speed, which makes it suitable for pa- 
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rameter search studies where a large number of models need to 
be run. It can also be coupled with the SMMOL radiative trans- 
fer code dRawlings & Yatesl 120011) for a detailed treatment of 
the PDR emission properties. Lee96mod was developed from 
the time-dependent chemical model by Lee, Herbst, and collab- 
orators. It is strongly geared toward complex chemical calcu- 
lations and only accounts for constant temperatures, neglecting 
local cooling and heating. Meijerink is a relatively young 
model with special emphasis on XDRs (X-ray dominated re- 
gions) which quickly evolved in th e course of this study and we 
refer to lMeiierink & Spaans (2005j) for a detailed review of the 
current status. In the Appendix we give a tabular overview of all 
main model characteristics. 



6. Concluding remarks 

We present the latest result in a community wide compara- 
tive study among PDR model codes. PDR models are avail- 
able for almost 30 years now and are established as a common 
and trusted tool for the interpretation of observational data. The 
PDR model experts and code-developers have long recognized 
that the existing codes may deviate significantly in their results, 
so that observers must not blindly use the output from one of 
the codes to interpret line observations. The PDR-benchmarking 
workshop was a first attempt to solve this problem by separat- 
ing numerical and conceptional differences in the codes, and re- 
moving ordinary bugs so that the PDR codes finally turn into a 
reliable tool for the interpretation of observational data. 

Due to their complex nature it is not always straightforward 
to compare results from different PDR models with each other. 
Given the large number of input paramters, it is usually possible 
to derive more than one set of physical parameters by compar- 
ing observations with model predictions, especially when one is 
chiefly interested in mean densities and temperatures. Our goal 
was to understand the mutual differences in the different model 
results and to work toward a better understanding of the key pro- 
cesses involved in PDR modeling. The comparison has revealed 
the importance of an accurate treatment of various processes, 
which require further studies. 

The workshop and the following benchmarking activities 
were a success regardless of many open issues. The major re- 
sults of this study are: 

- The collected results from all participating models rep- 
resent an excellent reference for all present PDR codes 
and for those to be developed in the future. For 
the first time such a reference is easily available not 
only in graphical form but also as raw data. (URL: 
|http : / /www . phi . uni-koeln . de/pdr-comparison| l 

- We present an overview of the common PDR model codes 
and summarize their properties and field of application 

- As a natural result all participating PDR codes are now better 
debugged, much better understood, and many differences be- 
tween the results from different groups are now much clearer 
resulting in good guidance for further improvements. 

- Many critical parameters, model properties and physical pro- 
cesses have been identified or better understood in the course 
of this study. 

- We were able to increase the agreement in model prediction 
for all benchmark models. Uncertainties still remain, visi- 
ble e.g. in the deviating temperature profiles of model V2 
(Fig. n~2b or the large differences for the H2 photo-rates and 
density profiles in model V4 (cf. online data archive). 



- All PDR models are heavily dependent on the chemistry and 
micro-physics involved in PDRs. Consequently the results 
from PDR models are only as reliable as the description of 
the microphysics (rate coefficients, etc.) they are based on. 

One of the lessons from this study is that observers should not 
take the PDR results too literally to constrain, for example, phys- 
ical parameters like density and radiation field in the region 
they observe. The current benchmarking shows that all trends 
are consistent between codes but that there remain differences 
in absolute values of observables. Moreover it is not possible 
to simply infer how detailed differences in density or tempera- 
ture translate into differences in observables. They are the result 
of a complex, nonlinear interplay between density, temperature, 
and radiative transfer. We want to emphasize again, that all par- 
ticipating PDR codes are much 'smarter' than required during 
the benchmark. Many sophisticated model features have been 
switched off in order to provide comparable results. Our inten- 
tion was technical not physical. The presented results are not 
meant to model any real astronomical object and should not be 
applied as such to any such analysis. The current benchmarking 
results are not meant as our recommended or best values, but 
simply as a comparison test. During this study we demonstrated, 
that an increasing level of standardization results in a significant 
reduction of the model dependent scatter in PDR model predic- 
tions. It is encouraging to note the overall agreement in model 
results. On the other hand it is important to understand that small 
changes may make a big difference. We were able to identify a 
number of these key points, e.g. the influence of excited hydro- 
gen, or the importance of secondary photons induced by cosmic 
rays. 

Future work should focus on the energy balance problem, 
clearly evident from the sometimes significant scatter in the re- 
sults for the non-isothermal models V1-V4. The heating by pho- 
toelectric emission is closely related to the electron density and 
to the detailed description of grain charges, grain surface recom- 
binations and photoelectric yield. The high temperature regime 
also requires an enlarged set of cooling processes. Another im- 
portant consideration to be adressed, especially when it comes 
to comparisons with observations is the model density structure, 
i.e. clumping or gradients. As a consequence we plan to con- 
tinue our benchmark effort in the future. This should include a 
calibration on real observational findings as well. 
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Appendix A: Characteristics of Participating PDR 
Codes 

In Tab. I A. 1 1 we summarize the most important characteris- 
tics of the participating PDR codes. This table summarizes 
the full capabilities of the PDR codes and is not limited 
to the benchmark standards. It has been extracted from de- 
tailed characteristics sheets, available online for all codes: 
|http : //www .phi . uni-koeln . de/pdr-comparison| 
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Table A.l. Full capabilities of the PDR model codes participating in the Leiden comparison study 
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photoionization 


X 


X 


X 


X 








X 






X 


carbon ionization heating 


X 


X 


X 


X 








X 




X 


X 


other species (Si, etc.) 


X 




X 


















gas-grain collisions 


X 


X 


X 


X 








X 




X 


X 


turbulence heating 


X 




X 


X 
















chemical balance 






X 


X 








X 









UV TRANSFER 
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Cloudy 


COSTAR 


Meudon 


Pi 

Q 

1 

u 


HTBKW 


■ 

< 

«5 

o 


Aikawa 


Leiden 


Lee96mod 


Sternberg 


Meijerink 


solved self-consistently 


X 


X 


X 


X 


X 


X 




X 




X 


X 


simple exponential attenuation 


X 


X 


X 


X 


X 


X 


X 


X 


X 


X 


X 


bi-exponential attenuation 




X 
















X 




full RT in lines 






X 










X 








DUST 


treatment of rad. transfer 


X 




X 




X 


X 




X 




X 


X 


grain size distribution 


X 




X 


X 




X 












extinction/scattering law 


X 


X 


X 


X 


X 


X 


X 


X 


X 


X 




albedo 


X 




X 


X 








X 






X 


scattering law 


X 




X 










X 








H 2 SHIELDING 


shielding factors 


X 


X 






X 




X 




X 


X 


X 


single line 


X 






X 














X 


detailed solution 1 


X 




X 






X 




X 








CO SHIELDING 


shielding factors 


X 


X 




X 


X 


X 


X 


X 


X 


X 


X 


single line 


X 














X 






X 


detailed solution 






X 










X 








isotope selective photodissociation 






X 






X 




X 






X 


UV PROFILE FUNCTION 


Gaussian 








X 


X 














Voigt 


X 




X 










X 




X 


X 


Box 
























other 
























RADIATIVE TRANSFER IN COOLING LINES 


escape probability 


X 


X 


X 


X 


X 


X 


X 


X 




X 


X 


other 
























IR pumping 


X 


X 






X 






X 






X 


OBSERVATIONAL LINES 


self-consistent treatment with cooling 


X 






X 
















escape probability 


X 






X 




X 


X 


X 




X 


X 


other 






X 






X 












H 2 


X 




X 










X 




X 




HD 






X 










X 




X 




ll CO 


X 




X 


X 


X 


X 


X 








X 


li co 


X 




X 






X 










X 


c l8 o 






X 






X 












U C 18 






X 






X 












[OI] 


X 




X 


X 


X 


X 


X 


X 




X 


X 


[CH] 


X 




X 


X 


X 


X 


X 


X 




X 


X 


[CI] 


X 




X 


X 


X 


X 


X 


X 




X 


X 


Si + 


X 




X 




X 












X 


CS 






X 
















X 


H z O 










X 














h^o 
























HCO + 






X 




X 


X 










X 


OH 










X 














[Sil] 


X 








X 


X 












[sinsii] 


X 








X 


X 










X 
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>, 

T3 


TAR 


s 



T3 


Pi 



°-, 


KW 


MA-t 


e8 


s 

Oi 


6mod 


iberg 


erink 




3 


t/2 


3 




0a 






3 


©\ 


C 






O 





t> 


u 


H 







"33 


tu 


a> 


at 




U 


u 


s 




S3 






J 




+j 

cc 




[Fel], [Fell] 


X 








X 












X 


COMPUTED LINE PROPERTIES 


resolved line profile 






X 




X 


X 




X 








continuum rad./rad transfer in UV 


X 




X 


















line center intensities 


X 




X 






X 




X 








line integrated intensities 


X 




X 


X 


X 


X 




X 






X 


optical depths 


X 




X 


X 


X 


X 




X 






X 


Gaussian line profile 


X 




X 




X 


X 




X 






X 


box line profile 
























turbulence included 


X 




X 






X 










X 


COLLISIONS 


H-H 


X 








X 










X 




H 2 -H 


X 




X 


X 


X 






X 




X 


X 


H 2 -H+ 


X 




X 














X 




H 2 -e 


X 




X 




X 










X 




H2 - H2 


X 




X 




X 






X 




X 


X 


CO-H 


X 


X 


X 


X 


X 






X 






X 


CO-H 2 


X 


X 


X 


X 


X 


X 




X 






X 


CO-e 


X 


X 




X 




X 










X 


CO -He 






X 


X 














X 


C-H 


X 


X 


X 


X 


X 


X 




X 




X 


X 


C-H 2 


X 


X 


X 


X 


X 


X 




X 






X 


C-e 


X 




X 


X 














X 


C-He 


X 




X 


X 
















C-H 2 
























C + -H 


X 


X 


X 


X 


X 






X 






X 


C + -H 2 


X 


X 


X 


X 


X 


X 




X 






X 


C + -e 


X 


X 




X 


X 


X 




X 






X 


O-H 


X 


X 


X 


X 


X 


X 




X 




X 


X 


0-H 2 


X 


X 


X 


X 


X 


X 




X 




X 


X 


-H+ 


X 




X 


X 








X 








O-e 


X 


X 


X 


X 


X 














O-He 


X 




X 


X 
















OH-H 
























OH -He 
























OH-H 2 










X 


X 












H-H 


X 






















H 2 0-e 
























H 2 0-H 






















X 


H 2 - H 2 










X 












X 


H 2 0-0 
























dust - H/H 2 


X 








X 














dust-any 


X 






















Si + - H 


X 




X 






X 












HD-H 






X 


















HD - H 2 






X 


















PAH-any 


X 








X 














OUTPUT 


abundance profile over (Ay/depth) 


X 


X 


X 


X 


X 


X 


X 


X 


X 


X 


X 


column density over (Ay/depth) 


X 


X 


X 


X 




X 








X 


X 


temperature profile over (Av/depth) 


X 


X 


X 


X 


X 


X 


X 


X 




X 


X 


emitted intensities 


X 




X 


X 


X 


X 




X 




X 


X 
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TAR 


a 

■d 


Pi 


1 




MA-t 


> 


s 


6mod 


iberg 


erink 




3 




3 




pa 






IS 


©\ 


C 









O 




u 


H 


O 




"33 






at 




U 


U 


s 




a 




<i 


-1 




+j 

cc 


£ 


opacities at line center 






X 




X 


X 




X 




X 


X 


heating and cooling rates over (Ay/depth) 


X 




X 


X 




X 




X 




X 


X 


chemical rates over (Ay/depth) 






X 


X 




X 




X 




X 


X 


excitation diagram of H2 


X 




X 










X 




X 





